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We calculate inclusive hadrons production in pA collisions in tiie small-x saturation formalism 
at one-loop order. The differential cross section is written into a factorization form in the coordi- 
nate space at the next-to-leading order, while the naive form of the convolution in the transverse 
momentum space does not hold. The rapidity divergence with small-x dipolc gluon distribution of 
the nucleus is factorized into the energy evolution of the dipole gluon distribution function, which 
is known as the Balitsky-Kovchegov equation. Furthermore, the coUinear divergences associated 
with the incoming parton distribution of the nucleon and the outgoing fragmentation function of 

^^", the final state hadron are factorized into the splittings of the associated parton distribution and 

I ■ fragmentation functions, which allows us to reproduce the well-known DGLAP equation. The hard 

^— V ' coefficient function, which is finite and free of divergence of any kind, is evaluated at one-loop order. 
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I. INTRODUCTION 



^v Inclusive hadron production in pA collisions have attracted much of theoretical interests in recent years |1|-|13||. 

f*\| In particular, the suppression of hadron production in the forward dAu scattering at RHIC observed in the experi- 

ments [3,[ia| has been regarded as one of the evidences for the gluon saturation at small-x in a large nucleus |3,|8„ 16]. 
Saturation phenomenon at small- a; in nucleon and nucleus plays an important role in high energy hadronic scatter- 
p^ ' ing [l7H20J. In this paper, as an important step toward a complete description of hadron production in pA collisions 
I , in the saturation formalism, we calculate the one-loop perturbative corrections. Previous attempts have been made in 
Mh the literature. In particular, in Ref. 5], part of one-loop diagrams were evaluated. However, the rapidity divergence 
is not identified and the collinear evolution effects are not complete. Recently, some of the higher order corrections 
were discussed in Ref. [9[, where it was referred as "inelastic" contribution. In the following, we will calculate the 
complete next-to-leading order (NLO) corrections to this process in the saturation formalism. A brief summary of 
\^ . our results has been published earlier in Ref. [21j|. 
Inclusive hadron production in pA collisions. 
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• , can be viewed as a process where a parton from the nucleon (with momentum p) scatters on the nucleus target 

^n . (with momentum Pa), and fragments into final state hadron with momentum Ph- In the dense medium of the large 
nucleus and at small- a;, the multiple interactions become important, and we need to perform the relevant resummation 
to make the reliable predictions. This is particularly important because the final state parton is a colored object. 
Its interactions with the nucleus target before it fragments into the hadron is crucial to understand the nuclear 
effects in this process. In our calculations, we follow the high energy factorization, also called color-dipole or color- 
glass-condensate (CGC), formalism [20, [22, [231 to evaluate the above process up to one- loop order. We notice that 

y^ , alternative approaches have been proposed in the literature |10l - [l3j to calculate the nuclear effects in this process. 

^_' According to our calculations, the QCD factorization formalism for the above process reads as, 
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-^ — ^xfa{x,^i)Dh/ciz,ii) / [dXA_]S'^,.{[xA)'Ha^cias,^,[xA_]^J,) , (2) 



where ^ = t/xz with r — p^e^ / ^, y and p±^ the rapidity and transverse momentum for the final state hadron and 
s the total center of mass energy square s = {p + Pa)"^, respectively. Schematically, this factorization is illustrated in 
Fig. 1, where the incoming parton described by the parton distribution fa{x) scatters off the nuclear target represented 
by multiple-point correlation function ^^([xj^]), and fragments into the final state hadron defined by the fragmentation 
function D}^i^{z). All these quantities have clear operator definitions in QCD. In particular, fa{x) and Di^/^{z) are 
collinear parton distribution and fragmentation functions which only depend on the longitudinal momentum fraction 
X of the nucleon carried by the parton a, and the momentum fraction z of parton c carried by the final state hadron h, 
respectively. From the nucleus side, it is the multi-point correlation functions denoted as S^^{x±) (see the definitions 
below) that enters in the factorization formula, depending on the flavor of the incoming and outgoing partons and 
the gluon rapidity Y associated with the nucleus: Y « lii{l/xg) with Xg being longitudinal momentum fraction. 

At the leading order, S^ represent the two-point functions, including the dipole gluon distribution functions in the 
elementary and adjoint representations for the quark and gluon initialed subprocesses |20|, respectively. Higher order 
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FIG. 1. Schematic plot of the factorization, where T-i indicates the hard factor, TZ represents the rapidity divergence which 
are factorized into the dipole gluon distribution of the target nucleus {A), Cp and C/ stand for the coUinear divergences which 
are absorbed into the parton distribution functions of the projectile proton (P) and hadron (Ph) fragmentation functions, 
respectively. 

corrections will have terms that depend on the correlation functions beyond the simple two-point functions. Because 
of this reason, the integral [dx±] represents all possible integrals at the particular order. 

To evaluate NLO corrections, we will calculate the gluon radiation contributions. At one-loop order, the gluon 
radiation will introduce various divergences. The factorization formula in Eq. ([5]) is to factorize these divergences into 
the relevant factors. For example, there will be collinear divergences associated with the incoming parton distribution 
and final state fragmentation functions. In addition, there is the rapidity divergence associated with ^^([xj,]). These 
divergences naturally show up in higher order calculations. The idea of the factorization is to demonstrate that these 
divergences can be absorbed into the various factors in the factorization formula. After subtracting these divergences, 
we will obtain the hard factors T-La^c, which describes the partonic scattering amplitude of parton a into a parton c in 
the dense medium. This hard factor includes all order perturbative corrections, and can be calculated order by order. 
Although there is no simple fc_L-factorization form beyond leading order formalism, we will find that in the coordinate 
space, the cross section can be written into a nice factorization form as Eq. ^. Besides the explicit dependence on 
the variables shown in Eq. ([5]), there are implicit dependences on pj_[x_l] in the hard coefficients as well. 

Two important variables are introduced to separate different factorizations for the physics involved in this process: 
the collinear factorization scale n and the energy evolution rapidity dependence Y. The physics associated with fx 
follows the normal collinear QCD factorization, whereas the rapidity factorization Y takes into account the small- 2: 
factorization. The evolution respect to u is controlled by the usual DGLAP evolution, whereas that for S^ by the 
Balitsky-Kovchegov (BK) evolution [23[[2J]- In general, the energy evolution of any correlation functions can be given 
by the JIMWLK equation[23|, and the resulting equation is equivalent to the BK equation for dipole amplitudes. 
In particular, our one-loop calculations will demonstrate the important contribution from this rapidity divergence. 
Schematically, this factorization is shown in Fig. [T] 

Our calculations should be compared to the similar calculations at next-to-leading order for the DIS structure func- 
tions in the saturation formalism [2Ch:2S| . All these calculations are important steps to demonstrate the factorization 
for general hard processes in the small-x saturation formalism [29]. The rest of the paper is organized as follows. In 
Sec. II, we discuss the leading order results for inclusive hadron production in pA collision, where we also set up the 
framework for the NLO calculations. Sec. III. is divided into four subsections, in which we calculate the NLO cross 
section for the q — >■ 5, g — > g, q — >■ g and g — > q channels . The summary and further discussions are given in Sec. IV. 



II. THE LEADING ORDER SINGLE INCLUSIVE CROSS SECTION. 



The leading order result was first formulated in Ref. jl|. For the purpose of completeness, we briefly derive the 
leading order cross section to set up the baseline for the NLO calculation. Let us begin with the quark channel in pA 
collisions. As illustrated in Fig. [51 the multiple scattering between the quark from the proton and the dense gluons 




FIG. 2. Typical Feynman diagrams for the leading order quark production qA 
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inside the nucleus target can be cast into the Wilson line 
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with A'^{x^ , x±) being the gluon field solution of the classical Yang-Mills equation inside the large nucleus target. 

Therefore, the leading-order cross section for producing a quark with finite transverse momentum k± at rapidity y 
in the channel qA — > qX can be written as: 
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with x = ^e^ and x„ = -^e ^. The notation 



)y indicates the CGC average of the color charges over the nuclear 
wave function where Y ~ \iil/xg and Xg is the smallest longitudinal momentum fraction of the probed gluons, and 
is determined by the kinematics Q. Normally, we first compute the correlator {TrU{x±)U^yj_)) in the McLerran- 
Venugopalan model[i9| as the initial condition, and then we perform the energy evolution for the correlator which 
introduces the rapidity (Y) dependence. The energy evolution equation at small- a; for dense nucleus targets is the BK 
equation as we shall demonstrate later when we remove the rapidity divergence. When multiplied by the fragmentation 
function, the above result will lead to the differential cross section for hadron production in pA collisions. 

It is straightforward to include the gluon initiated channel, and the full leading order hadron production cross 
section can be written as 
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d^p^dyh 
with px — zk±, Xp = -^^e^'*, r = zxp and Xg — -^y=e^^''. Here we have defined 

Hk^)-l^^-^^^e-^'^^<^^-y^^S^^\x^,y^), 

with Sy {x±,y±) ~ — (TiU{x±)U^ {y±))y. T{k±) is defined similarly but in the adjoint representation 

Hk.)^ j'^^j0^e'^-^-i^^'y^^~S^\x.^y.), 

where 5y (xj_,yj_) = ^/_-^ (TrVF(x^)VF^(y_L)).j, and W{x) is a Wilson line in the adjoint representation. It represents 
the multiple interaction between the final state gluon and the nucleus target. In general, the adjoint Wilson lines can 
be replaced by two fundamental Wilson lines by using the identity 

W''^{xi_) ^ 2Tr [T''U{xi_)T^U\xi_)\ , (8) 
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^ Here we are only interested in tlie inelastic production of the quark in the forward scattering which produces quark with finite 
transverse momentum. There is also elastic scattering contribution to the cross section which generates vanishing fcx, such as 
y^, xqf{x)S^-^' {kj_) J d?b to the total cross section. 



and the color matrices can be removed using the Fierz identity T^jT^i — ^SuSjk — j^SijSki- It is straightforward to 
show that 
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i^^uyi-) = j^^ [{TrU{xi_)UHyi_)TrU{yi_)U\xi_))y - l] , (9) 



which, in the large Nc limit, allows us to write 

72^ . W2 



It is very important to keep in mind that the normalization of the dipole amplitudes S^'^'{x±,y±) is unity when 
x± = y±. In addition, since normally (TrU{x±)U^{y±)')y is real, it is easy to see that S'-^''{x±,y±) = S'-^\y±,x±). 



K{^±-y^r 



the 



If we further neglect the impact parameter dependence, one will find that S^'^'{x±,y±) = exp . 

McLerran-Venugopalan model, where Qs is the saturation momentum which characterizes the density of the target 
nucleus. The analytical form of the dipole amplitude can help us to test the properties of dipole amplitudes mentioned 
above. 

We would like to emphasize that in Eq. ([5]) we do not include the transverse momentum dependence in the incoming 
parton distribution from the nucleon. In the forward pA collisions, the transverse momentum dependence from the 
incoming parton distribution of the nucleon is not as important as that from the nucleus target. Therefore, in the 
current calculations, we neglect this effect. As a consistent check, the one-loop calculations in the following support 
this assumption. In particular, the coUinear divergence associated with the incoming parton distribution contains no 
transverse momentum dependence. 

III. THE NEXT-TO-LEADING ORDER CROSS SECTION 

In this section, we will present the detailed calculations for the NLO corrections to the leading order result in 
Eq. ([5]). There are four partonic channels: q — ^ qg, g — > gg, q — )■ gq, g — > qq. We will carry out the calculations for 
these channels separately. 

A. The quark channel q ^ q 

The quark production contribution contains the real and virtual gluon radiation at the NLO. For the real contribu- 
tion, we will calculate q ^ qg first. The real diagrams with a quark (with transverse coordinate b\) and gluon (with 
transverse coordinate xj_) in the final state, as shown in Fig. [3l have been studied in Ref. [30l - l32l |. We take eq.(78) 
of Ref. [34I as our starting point which giveqj 
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SP {b^,x^,b'^,x'J + S^y'> {vi_ , v'^ ) 
-S^' {b± ,x±,v'^) - S^' {v± , x'^ , b'^) 
where u± = x±_ — b±, u'j_ = x'^ — 6'^, v±_ = (1 — ^)x± + ^b±, v'j_ = (I — S,)x'^ + S,b'^ and 

S^^\b^,x^,b'^,x'^)^-^^(TT{U{b^)UHb'^)T''T^)[Wix^)WHx'^)Y')^, (12) 

S'^\b^,x^,v'J = ^ (Tr {Uib^)T^UHv'jT^) W^''ix^))y . (13) 



For convention reasons, we have interchanged the definition of z and 1 — z and replaced the variable 2 by ^. 
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FIG. 3. The real diagrams for the next-to-leading order quark production qA -^ q + X. 



For a right-moving massless quark, with initial longitudinal momentum p+ and no transverse momentum, the splitting 
wave function in transverse coordinate space is given by 
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{^a+5l3+ + i6a-5j3-), A = 2. 



(14) 



where A is the gluon polarization, a, p are helicities for the incoming and outgoing quarks, and 1 — ^ = ^ is 
the momentum fraction of the incoming quark carried by the gluon. Since the Wilson lines in the fundamental 
representation and the adjoint representation resum the multiple interactions of quarks and gluons with the nucleus 
target, respectively, one can easily see that these four terms in the last two lines of the Eq. (|lip correspond to those 
four graphs in Fig. [3] The Sy term which corresponds to Fig. |3] (a) and resums all the multiple interactions between 
the quark-gluon pair and the nucleus target, represents the case where interactions take place after the splitting both 
in the amplitude and in the conjugate amplitude. The Sy term which comes from Fig.[3](b), resums the interactions 

before the splitting only and the Sy terms represent the interference terms as shown in Fig.[3](c) and (d). 

There are two contributions for inclusive hadron production at the next-to-leading order, namely, quark productions 
associated with D^/q which is indicated by the cross in Fig. [3] (while the gluon is integrated) and gluon productions 
associated with the fragmentation function D^/g (while the quark is integrated) . 

Let us study the former case by integrating over the phase space of the final state gluon [k^ , fci^). We can cast 
the real contribution into 
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where x = t/z^ and Cp — {N^ — l)/2A^c, and I and J^ are defined as 
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FIG. 4. Typical virtual diagrams for the next-to-leading order quark production qA -^ q + X. 



'{Ti-[U{x±)U^b±)]TT[U{b±)U'' {y±)])Y ■ Several steps are necessary in deriving the above 
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and Sy {x±,b±,y±) = 

result from Eq. (jlip . By integrating over the gluon momentum, we identify x± to x^ which simplifies Sy to 
S^'^' {b±,b'j_). This is expected since we know the multiple interactions between the gluon and the nucleus target 
should cancel if the gluon is not observed. Furthermore, using the Fierz identity, one can write 



Sy ib±,X±,v'j^) 
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which only involves the Wilson lines in the fundamental representation. Then, the final steps, which include the Fourier 
transforms, as well as the convolutions of the quark distribution and fragmentation function, are quite straightforward. 

Before we proceed to the calculations of the virtual diagrams, we comment on the result shown in Eq. (|15p . The 
major obstacles of evaluating the integrals in Eq. dTSl) are the divergences. There are three types of singularities 
lying in that equation, namely, the rapidity divergence which occurs at ^ = 1 when the rapidity of the radiated 
gluon becomes — oo, and the coUinear singularities which correspond to the cases that the final state gluon is either 
collinear to the initial quark or final state quark. We shall expect that the virtual diagrams cancel some part of the 
divergences, while the uncancelled divergences shall be absorbed into the renormalization of the quark distribution 
and fragmentation functions as well as the target dipole gluon distribution {Sy {x±,y±)). After these subtractions, 
the remainder contributions should be finite and give us the NLO correction to the single inclusive hadron production 
cross section. 

The evaluation of the virtual graphs as shown in Fig. 2] are quite simple in the dipole picture. Their contributions 
are proportional to 
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where the factor of 2 takes care of the fact that the mirror diagrams of Fig. |4] give the identical contributions when 
the virtual loop is on the right side of the cut. It is straightforward to see that these two terms in the last line of 
Eq. P^ correspond to the Fig. 2] (a) and (b), respectively. This eventually leads to 
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where explicitly one writes 
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{q± ~ ^k^)Hq± - kgl±r 

It is easy to see that the virtual contributions indeed contain three types of singularities as we mentioned before. There 
are two important features that we wish to emphasize here. First, the rapidity divergence term is only proportional 
to Nc/2 since I vanishes at ^ — > 1 limit. This agrees with the BK equation since there is no 1/N^ corrections to 
the leading order BK equation. Second, when one integrates over the quark transverse momentum k±, the rapidity 
divergence disappears due to the complete cancellation between the real and virtual contributions. 



1. The rapidity divergence 

Now we are ready to evaluate NLO contributions by the following procedures. First, we remove the rapidity 
divergence terms from the real and virtual contributions by doing the following subtractions 
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where T^'^'{k^) is the bare dipole gluon distribution which appears in the leading order cross section as in Eq. ([5]) 
and it is divergent. T{k±) is the renormalized dipole gluon distribution and it is assumed to be finite. We can always 
decompose the dipole splitting kernel as 
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where the first two terms are removed from the virtual contribution while the last term is removed from the real 
diagrams. This procedure is similar to that for the collinear factorization, where we modify the bare leading order 
parton distributions to the finite parton distribution with the higher order radiation. Using Eqs. ([6| I21|) . we can see 

that the differential change of the dipole amplitude 5^ (a;_L, yj_) yields the BK equation 
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It is important to note that if we conduct the leading order classical calculation, we will not get any energy dependence, 
namely the Y dependence, in the scattering amplitudes. It is the BK evolution equation as shown above which gives 
the energy dependence to those scattering amplitudes. To derive the BK equation from Eqs. (|6l [2T1) . one needs to reset 



the upper limit of the d£, integral in Eq. ([2T|) to 1 — e , with Y being the total rapidity gap between the projectile 
proton and the target nucleus. Here y — > cx) as the center of mass energy s — )■ oo. By doing so, we introduce the 

(2) 

rapidity Y dependence, namely the energy dependence, of the two-point function Sy {x±,y±) from which the BK 
equation can be understood and therefore derived. Another way to derive this equation is to slightly move away from 
the light cone as in the derivation of the Balitsky equation[23|. The rapidity divergence is an artifact that we put 
both the projectile and targets on the light cone in the high energy limit. By slightly tilting away from the light cone, 
we can modify the ^ integral and obtain L ■_,j _y . In addition, when one integrates over the transverse momentum 
fcx as in Eq. (PT|) . one finds that the rapidity divergence disappears as expected (33| . 

The physical interpretation of the rapidity divergence subtraction is quite interesting. Although the soft gluon is 
emitted from the projectile proton which is moving on the forward light cone with the rapidity close to +oo, it is 
easy to see that the rapidity of this soft gluon goes to -co when ^ — >■ 1. As a matter of fact, this soft gluon can be 
regarded as collinear to the target nucleus which is moving on the backward light cone with the rapidity close to — cx). 
Therefore, it is quite natural to renormalize this soft gluon into the gluon distribution function of the target nucleus 
through the BK evolution equation. 

After the subtraction of the rapidity divergence, both of the real and virtual contributions become regulated in 
terms of the d^ integral which leads to the change of the splitting function into .j'"^! . Here we introduce the following 
property of the plus-function 
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where g{S,) can be any non-singular functions, while /(i^) is singular at ^ = 1 and (/(^)) , is regulated. 
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2. The collinear divergence 



The second step is to use the dimensional regularization {D — A — 2e) and follow the MS subtraction scheme, in 
order to compute and remove the collinear divergence from both real and virtual contributions. For convenience, we 



introduce the following integrals, 
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Clearly, there is no divergence in I^. Let us take the evaluation of /i(fc^) as an example. As standard procedure 



in the dimensional regularization (I? = 4 — 2e) and the MS subtraction scheme, we change the integral / ,^Jt^ into 
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M / (-9 -12^2^ where /^ is the scale dependence coming from the strong coupling g. Using Eq. ^ and the identity 
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where cq = 2e '*'^, 7^; is the Euler constant and r± — x± — y±. 
To evaluate l2{kj_), we first rewrite it as 
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where I21 is finite and defined as 
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The basic idea is to subtract a term which is proportional to ln(f — ^)~ from I2. It is quite straightforward to show 
that the last two terms in the above equation give J-{k±) ln(l — ^)2 by using the integral identity 
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To compute and remove the coUinear divergence in the virtual diagrams, one needs to use the following integral 
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where the usual Feynman integral trick is used in the derivation. Therefore, setting the quark distribution, the 
fragmentation function and the splitting function aside, the virtual contribution can be cast into 
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where l3v{ki_) is finite and defined as 
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{q±-k^y {q±-iki_Y{q^-kgiA_f 

^k^r 



2tt 



fci 



(33) 



To derive the above expressions, Eq. (1301) is used repeatedly. It is also useful to notice that 



d2fcij_e-*'=i-L-*^^ln 



(fci± - CkA 



A-K 



Sir±) 



Lkpi'si-'-l p-ii'k±-r± 



(34) 



which can lead us to the final factorized formula. 

By combining the collinear singularities from both real and virtual diagrams, we find the coefficient of the coUinear 
singularities becomes 'Pqq{S,) which is defined as 



^,g(C) = 



1-e 



1 + ^^ 3,,^ ^, 



(35) 



Now we are ready to remove the collinear singularities by redefining the quark distribution and the quark fragmentation 
function as follows 



q{x,^i)^q<^^\x) 



las{^i) fU^ 



e 2n 



^-CpV^qiOq,^ 



i?,/,(z,M) = <,(.) - \^f fc^v,,m,/, (0 , 



(36) 
(37) 



which is in agreement with the well-known DGLAP equation for the quark channel. We will be able to recover the 
full DGLAP equation once we finish the calculation on all channels. Using Eq. ([S]) and combine it with the NLO real 
and virtual contributions, it is almost trivial to show Eq. (|36)) . It is a little bit less trivial to derive Eq. ((37)) . By 
combining the relevant terms in the real and virtual contributions, we obtain a term which reads 



\'^ [^Dr,,,{z) j' ^d^C,Vqq{Oxq{x)^Tl^ 



t) 



By changing variable z' — z^, we can rewrite the above term as 

1 ctsi^j) f dz' 



e 2tt 



-,12 



xq{x)F (^) £ ^CpVq,{S,)D„,q [^-^ , 



(38) 



(39) 



which allows us to arrive at Eq. (|37|) easily by combining this term with Eq. ([5]). 

One might worry about the term which is proportional to -^F (fc^) {Cp — ^f^ ln(l — ^)^ since it is logarithmically 
divergent when ^ — )• 1 . Let us show that this singularity will also cancel between the real and virtual contributions as 
follows 



^^7i T^^^i^) ln(l - 0^ - Xpq{xp) / di- -— 



in(i - e 



;.(<i±i^lM!)^.«, 



(40) 



where the first term on the left hand side of the above equation comes from the real diagrams while the second term 
comes from the virtual graphs. Here we have used Eq. (j24l) again. 



3. Finite contributions 



Now we have removed all the collinear singularities by renormalizing the quark distribution and the quark fragmen- 
tation function. The rest of the contribution should be finite. The last procedure is to assemble all the finite terms 
into a factorized formula. For the quark channel contribution: qA -^ h -\- X^ we find that the factorization formula 
can be explicitly written as 



dyd'^p± 



= I — — ^xq{x,^)Dh/qiz,fi) 



d'^x±d^yA 



{2^y 



[s'^\xu 



yi 






iM 






(41) 



10 
up to one- loop order. The leading order results have been calculated as shown in Eq. ([5]), from which we have 

<i = e-'"-''"'^(l-0, (42) 

where k± = p±/z and r± — x± — y±. Our objective here is to compute the hard coefficients 'H2n„ ^^'^ T~^\qq- It is just 
straightforward to show that "HLi reads as follows 



<l = CFP,,(01n-|^ 



r^/i^ 



ik\ ■r\ 



iCpSil-Oe 



-ik± -rx 



In. 



r^ 
Cq 



{2Cf ~ Nc) e 



—ik± ■r± 



-J21 — 



(1-0- 



i-e 



(43) 



where the terms in the first line come from the finite logarithmic terms in Ii(k±) and Iy{k±), and /21 is calculated 
from l2i{k±) which yields 



l2l = 



d%. 



(^6i - 



bl i^b^ - r^) 



1 



„ — ik I -b I 



1 



(44) 



It is clear that the last term comes from the ln(l — ^)^ terms as we have shown in Eq. (|40l) . It is also important to 
note that the second line in Eq. (j43|) (the l2{k±) type term) drops out if we take large Nc limit. The large Nc will 
greatly simplify our calculation in many aspects as we will show in the following sections. Furthermore, by choosing 
H = co/r± for the factorization scale, we can further simplify the above expressions. In the end, only the last term in 
the first line of the Eq. (|43)) survives. Since r± is of the order 1/Qs in the saturation regime, one can easily see that 
the factorization scale fi c:i Qs in terms of the above choice. 

The second hard coefficient H^J^ is related to the non-linear terms such as /3(fc_L) and l3v{k±) which give 



0/(1) ^ 

''-41313 



-ik^-r^ ] ^^zi^k^.(x^-b^) l+C'^ 1 XjL-bA_ ^ yjL-b± 



-AirNre 



r^ 1 + £'^ 



'g-Hl-e)k±-iv±-b^) 



{bA 



yjL) 



-5^^\b^-y^) hfr', 



(45) 



where the first and second term in the curly brackets are calculated from l3{k±) and /3t,(fe±), respectively. 

To summarize the above results, we have demonstrated the QCD factorization for inclusive hadron production in 
the quark channel of pA collisions in the saturation formalism, and we have computed the NLO cross section in this 
processes. Clearly, the naive form of the k± factorization formula, which involves the convolution of unintegrated 
gluon distributions in the transverse momentum space, does not hold. Other channels can be calculated accordingly 
following the same procedure. 

4. The McLerran-Venugopalan model 

In this subsection, we calculate the hard coefficients in the well-known McLerran-Venugopalan (MV) model [13, H, 
|35| . In terms of the phcnomcnological application with additional parametrization of the saturation momentum, it 
is also known as Golcc-Bicrnat-Wusthoff (GBW) model [36]. In the MV and GBW model, if we neglect the impact 
parameter dependence for the sake of simplicity, the dipole scattering amplitude is parametrized as 



5'mv(2;-L,2/±) =exp 



ix^-yi^?Ql 



(46) 



which leads to J^{q±) 



ttQ\ 



■ exp 



. 9i 



where S±_ is the transverse area of the target hadron. In addition, if we 



further take the large Nc limit, the integral (Px±(Py±cPb± can be performed explicitly, which leads to the differential 
cross section depending on p± and Qs , 



dycPpi 



X 



— — £.xq{x,n)Dh/q{z,n) 



n 



(0) 
2qq 



"«a}(i) 



2tt 4w. 



(47) 



where 



n 



(0) 
2qq 



-^a-o^^exp 



•x/(i) ^^ 






In 






11 



(48) 



'^2 2 



^,«(0-^ 



n 



(1) 

4g9 



^ (i-e)+fci 

+7Vc5(l-0-^(fc±) 



In 



In 



+ exp 

_Ql_ 

//2e7E 
0^ , 






L(i'O) -1,- 



fc2 



--(A)^""^(-' 






fcj_eTE 



exp 
r)2 



h exp 






-(1,0) 



t: In , 9 — 



1 — exp 



-1,- 



fc2 






exp 



'2 1,2 



C'k 



L 



(1:0) 



/2l,2 



-1 



C'k 



(49) 



(50) 



where kj_ — pj_/z as above and L^^'^' (— 1, — x) = — [7^ + F (0, — a;) + log (— x)] e~^ is the Multivariate Laguerre 
Polynomial. L'^^'^' ( — 1,— x) is zero at a; = 0,oo, and reaches its maximum at x around 2. 

An important aspect of the above results is that we can compare with the coUinear factorization results in the dilute 
limit. For example, the forward quark production in pp collisions is dominated by the i-channel qg -^ qg subprocess in 
the coUinear factorization calculation. Because of the i-channel dominance, we find that the differential cross section 
can be written as 



(Pa{pp 



(Pk±dy 



for\ 



id limit 



—^IKXj-^^x g{x ) 



(51) 



in the forward limit, where fc_L and y are the transverse momentum and rapidity of the final state quark, respectively. 
The above result was obtained by taking the limit of — t ^ s ~ —u for the Mandelstam variables in the partonic cross 
section. Here, q{x) and g{x) are quark and gluon distributions from the incoming two nucleons, respectively. 

As a consistency check, we can take the dilute limit which gives fc^ ^ Q^, and obtain the leading contribution of 



Eq. (j47|) which reads 



^3^p+A-^h+X 



dycPpi 



dilute 



dz 



1-^ 



xq{x,n) 



a, 2N,S^Ql 



27r 



7rfc4 



(52) 



In arriving at the above result, we have also taken the limit ^ — >■ 1 (note that ^ 7^ 1 for real contributions due to 
subtraction) which corresponds to the limit — i ^ s ~ —u. We further notice that the quark saturation momentum [34l . 
1371 Q^g — "^^"° ^R^ — h'^px'G{x') with x'G{x') corresponding to the gluon distribution in a nucleon, p being the nuclear 
density, R being the size of the target nucleus and b being the impact parameter. In the dilute regime, the gluon 
distribution is additive in the target nucleus which allows us to write x'Ga{x') = 2\J B? — b'^S±px'G{x') = Ax'G{x') 
with A being the nuclear number |3 At the end of the day, by setting yzi = ^ which recovers the integration over 
the gluon longitudinal momentum fraction, we can obtain 



d^aP+^- 



■h+X 



dyd^PA 



dilute 



dz 

— Dh/q(z,p) 



dx 



-xq{x,n) 



2alx'GA{x') 



(53) 



which agrees with the coUinear factorization result for the quark channel. The comparison for all other channels shall 
follow in the same way. In conlusion, the factorization in Eq. ([2]) is consistent with the coUinear factorization result 
in the dilute limit in the forward pA collisions. 

B. The gluon channel g ^r g 

The computation for the g -^ g channel is very similar to the calculation we have done for the q ^ q channel. 
However, there is an additional complication in this calculation. As we will show later in the detailed derivation. 



Rigorously, one should write 5x = J d^b and use the relation that p f cPb2\/ K^ — b'^ 



p^R' 



12 




'ooooGeooooooo( 




{b) 



'U6m(WS66S6<^k^iQWSmi 




moGmcmGOGi^k^wsiosiQwssaQ., 




(c) 



id) 



FIG. 5. The real diagrams for the next-to-leading order gluon production gA -^ g + X. 

the sextupole, namely the correlation of six fundamental Wilson lines in a single trace, will start to appear in 
the cross section. The small- a; evolution equation of sextupoles [38| is different from the well-known BK equation 
which is derived for dipoles. This is normal since the quadrupoles also follow a different version of small- a; evolution 
equation[30|,|39|. The numerical study of the evolution for sextupoles is not yet available. Fortunately, the contribution 
from sextupoles is suppressed by a factor jj^ as compared to other terms. In addition, in principle, the four-point 

function S^'^^x±,b±,y±) can not be factorized into S^^^x±,b±)S'^'^'){b±,y±) unless the large Nc limit is taken. By 
taking the large Nc limit, not only can we simplify the calculation significantly, but also we can show that all the 
relevant S'-matrices are dipole amplitude S*"' which is universal at both leading order and NLO. From the universality 
point of view, it seems that the large Nc limit is essential to the factorization. Therefore, in our following derivation, 
we will take the large Nc limit right away, but we will comment on the property of the Nc corrections. 

The real diagrams, as shown in Fig.[5l have been studied in Rcf. [32]. Let us first analyze the S'-matrices associated 
with each graph in Fig.[5l For Fig.[5ja), before we integrate out the phase space of the unobserved gluon, we find that 

the multiple scacttering gives / fade [W{x±)W^ {^'±)] [W{b±)W^ (^1)] ^'^ fate ) , where x± and x'j_ are the transverse 

coordinates of the observed gluon in the amplitude and complex conjugate amplitude, respectively. Here, b± and 6^ are 
the coordinates of the unobserved gluon. By integrating over the phase space of the unobserved gluon, we identify b± to 
b'j^ which allows us to greatly simplify the above expression and obtain Nc [(TiW{x±)U{x'^)TrW(x'j_)U{x±)')y — l). 
The interaction between the unobserved gluon and the nucleus target is cancelled as expected. By taking the large 
Nc limit, we can further drop the second term and factorize the results into NcS^'^^x±,x'jJS^^''{x'j_,x±), where a 
factor of -^ has been attached as the color averageQ. Similarly, the Fig.[5{b) yields NcS^^'>{v±,v'^)S^'^\v'j^,vj_) with 
£,x± + (1 — £,)b± and v'j_ = S,x'^ + (1 — S,)b± Q For Fig. Eljc), we find the scattering matrix is proportional to 



v± 



{fadeW'"'{x^)W^''{b^)fftcW''f{v'^))y 

{TrUHv'^)U{xi_)TTUHx^)U{b^)TrUHb^)U{v'^))y 
- {TYUHx^)U{v'^)U^{b^)U{x^)U\v'^)U{b^))y , 



(54) 



where we have used Eq. ([5]) and if'^'"^T'' — [T", T''] in the derivation. In addition, we have assumed that the expectation 
value of the Wilson lines is real, which allows us to get, for example. 



{TrUHv'^)U{x^)TrU^ix^)U{b^)TrU^ib^)U{v'^))y 
- {TYU\x^)U{v'^)TYU^{v'jU{b^)TTU\b^)U{x^))^ . 



(55) 



* Strictly speaking, this factor should be ^_ since the number of gluon color is N^ — 1. hi the large Nc limit, we just put it as A^^. 

^ The way that we choose to define v_i_ and d^ here is to put the rapidity divergence at ^ = 1 according to the convention that the 
unobserved gluon's longitudinal momentum becomes infinitely soft. 
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FIG. 6. Typical virtual gluon loop diagrams for the next-to- leading order gluon production gA -^ g + X. 

The last term on the right hand side of Eq. ([54| is the sextupole that we discussed earlier and it is suppressed by -^ 

as compared to the first term. It is easy to see that the first term is proportional to N^ since it has three color traces. 
Therefore, we obtain that Fig. [Sljc) gives NcS^^^xj_,v'j_^)S^'^^v'j^,b±)S^^^b^,x±) in the large Nc limit. Similarly, 
following the same procedure, we find that Fig. [SJd) yields NcS^'^^v±,x'j_)S^^\x'^,b±)S^'^'>{b±,vj_). 

Now we can follow Ref. [32| and write down the cross section of producing a hadron with p± at rapidity y from a 
gluon as follows 



, pA^hX 



Aa/3 

X 's^^\x^,x'^)S'^^\x'^,x^) + S^^\v^,v'^^S'^^\v'^,vs_) 
~S^^\x^,v'^)S^^\v'^,b^)S^^\b^,x^) 
-S^^\v^,x'^)S^^\x'^,b^)S^^\b^,v^ 

where the g ^ gg splitting kernel is found to be |j 



E<*"/3(?,"lXc./3(^,"±) = 4(27r) 

\al3 



1-? 



+ e(l-0 



7? ? ■ 



(56) 



(57) 



with ux = x±_ — b± and u'^ 

under the interchange ^ O (1 — ^), and can be simply written as * ~5i ~^V ■ It is clear that the real contributions 

contain the rapidity divergence at ^ — > 1 limit. 

Similar to the quark channel, the virtual gluon diagrams as shown in Fig. [S] can be calculated accordingly, and we 
obtain 



x'j_ —b±. In addition, we find that the ^ dependence of the splitting function is symmetric 



^ dz f \ [^ jt f d^v±_ (fv'^ (fu± 

^D,/,{z)xM^p) I d^J (2^)2 (2^)2 (2^)2 



xe 



—ik± ■{v± — 



"^^ E 'f'a^cM^ ^' uj_)i:^^^p{p+,t u^) 



Xap 



' S(^\v^,v'^)S(^\v'^,v^)- S(^\b^,x^)S^^\x^,v'JS(^\v'^,b^) 



(58) 
(59) 

(60) 



where the factor of 2 comes from the fact that the mirror diagrams give the identical contributions when the virtual 
loop on the right side of the cut, while the factor of ^ is the symmetry factor arising from two identical gluons in the 
closed gluon loop. The virtual contribution contains rapidity divergence when ^ approaches and 1. This is easy to 
understand since the virtual gluon loop contribution is symmetric under the interchange ^ f^- (1 — ^). Assuming that 
S^'^i{x±, x'^) = S^'^'{x'^, x±) and use x± = v± + {l — ^)u± and b± = v±—^u±, one can easily show that the last line of is 



Here we have included the factor of -rp, which is in the spUtting kernel, into the cross section. 
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FIG. 7. Typical virtual quark loop diagrams for the next-to-leading order gluon production gA -^ g + X. 



symmetric under the interchange i^ O (1 — ^). Therefore, we can rewrite the splitting function j^ + ij^ + ^(1 — ^) 

as 2 Y3T -|- 2^(1 ~ fo^' the virtual part. Now the virtual contribution only contains rapidity singularity at ^ = 1. 

Following the procedure that we have illustrated above for the quark channel, we remove the rapidity divergence 
terms from the real and virtual contributions by doing the following subtractions 



jr(A:^)=jr(o)(fcj 



a,Nr. 






S'^'\x^,yj_)S'^'\y, 



J-,X_l] 



d'^x±d'^y±(Pb± 



-ik±-(x±-y±) _ 



ix± - yA_Y 



{x^-h^f{y^-h^Y 



s^^Hx±,y±)s'^^Hy±,b±)s'-^Hb 



-LiX±, 



(61) 



where T^'^' {k± ) is the bare dipole gluon distribution in the adjoint representation which appears in the leading order 
cross section as in Eq. ([SJ and it is divergent. J-{k±^ is the renormalized dipole gluon in the adjoint representation 
distribution and it is assumed to be finite. To arrive at Eq. (j6ip . we have taken the large Nc limit which allows us to 
neglect the sextupole and constant term which are suppressed by -^ . the full subtraction should include those terms 
as well. 

Now we are ready to show that Eq. (|5T|) is equivalent to the adjoint representation of the BK equation. The 
non-linear small-x evolution equation for a color dipole in some arbitrary representation R can be found in Eq. (5.18) 
in Ref. [40|. This equation reads 






(triiKW), 



d'^z± (x_L - j/±)^ 

{xi_ - Zj_)2(yj_ - 2_l)^ 



(62) 



where V is the Wilson line in the i?-representation. If one takes the fundamental representation, one can easily recover 
the BK equation as shown in Eq. ((23)) . If one sets V = W and uses the adjoint representation for the color matrices 
f^^ = —ifabc, one can use Eq. ([S]) to convert everything into the fundamental representation. It is straightforward to 
find Cr — Nc and 



{tTAW^^Wy^)^ = (TrC/t(a;^)f/(y^)Trt/t(y^)C/(a;^)), 



1 



- {TYU\x^)U{y^)U\z^)U{x^)U^{y^)U{z^))y , 



(63) 



where we have also assumed all the correlation functions on the right hand side of the above equation are real. By 
putting above expressions into Eq. (1621) . we can obtain the adjoint representation of the BK equation which is in 
complete agreement with Eq. ([CTjl if one also includes the large Nc corrections in Eq. (pTj) . This version of the BK 
equation actually contains the sextupole correlation term and constant term which coincide with Eq. (|54p and the 
discussion above. One can see that the cancellation of the rapidity divergence is complete, even if one includes all the 
large Nc corrections. After the subtraction of the rapidity divergence, the splitting functions become regulated and 
we can replace j^ by /^K) ■ 

Furthermore, before we take care of the coUinear singularities, we should also compute the quark loop virtual 
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diagrams as shown in Fig. [71 



9 Airr /' ^^n ( \ ( ^ [\c [ '^'^^ '^'^ '^'"^ 
-2a,N,T^, / -,Du„{z)x,9{x,) j^ d^ j (,^), (,^), (,^), 



xe 



— 'i/cj_ -{v^ - 



""^ E <*aMP+, ?, Ui)V^,-„;3b+, e, "±) 



Aq/3 



5(2)(„^, ^;;j_)5(2)(^^^ „^) _ s^^Kx^,v'^)S^^\v'^M) 



where T^ = g ^i^d the g —^ qq sphtting kernel is found to be 



E<l«/^(P^'^'"^X-"/3(P'''^'"^) = 2(27r)2 [e^ + (1 - 0'] 4- 



(64) 



(65) 



Aq^ 



Normally the g ^ qq channel is suppressed by one factor of l/N^ as compared to other leading N^ channels. However, 
the quark loop gains a factor of Nf since different flavors of quarks can enter the virtual quark loop. Nf is usually 
taken to be 3 which is the same as N^- Therefore, we also compute this channel since this might be as important as 
other contributions in terms of the numerical studies. There is no rapidity divergence in the quark loop contribution, 
however, it does contain collinear divergence. 

To compute and remove the collinear divergence, we define 



JC{k±,l±,q±) = j 



d^x±d^b±d y±d X' 



2^1 



(27r)6 



A^ -ik±-{x±-b±)-il±-(b±-y±)-iq±-(y±-x'^) 



(66) 



where the variable x'j^ is redundant and can be integrated out to give the area of the target nucleus, if one neglects 
the impact factor dependence. This allows us to transform Eq. (|56p into 



dal„7'^^ a.Nr f^ dz 



real 



d^p±dy TT^ 



:2^h/giz) 



.[1-^(1-0]^ 



;/C(gij_,g2J.,g3J 



k± - qi± + q3± 



xg{x) / d^qi^d^q2±d^qs.i_ 



k± - ^qi± + £,q2± 



(fc± -gi± + 93±)^ (fc± -^gi± +^92±)^ 



' (fcj_-IJij_+(j3_ 



It is quite clear that among those three terms in the above equation 

which should be absorbed by the gluon distribution, j-, j ^ r^ 

should be associated to the fragmentation function, while the crossing term is finite. Similarly, the virtual gluon loop 
contribution can be transformed into 



w term gives the collinear singularity 
w term yields the collinear singularity which 



dz f^ 

—Dh/g{z)xpg{xp) I d£, 

z Jo 



(T^"5«'-« 



X / d'^qi±d'^q2j_JC{qi_i_,q2A_,q2± - k±_) / d'^l± 



l± - S.k± + q2± - gi. 



ll il±-^k^ + q2±-qi±y 



(67) 



and the quark loop term can be turned into 

27r2 



D,,^g{z)xpg{xp) / d^ [e + (1 - e)'] 



d'^qi±G{qi±,qi± - ki_) / d'^li_ 



?fcj 



9iJ- 



ll {l±-£,k±+qi±y' 



(68) 



It is straightforward then to use Eq. pip to compute collinear singularity for the virtual contributions. 

By combining the collinear singularities from both real and virtual diagrams, we find the coefficient of the collinear 
singularities becomes Vgg{S,) which is defined as 



VggiO 



^ 



1-e 



+ e(i-c) 



11 2NfTR 
~6 3Nc 



-5(1-0, 



(69) 
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where the first term comes from the real diagrams, the term which is proportional to ■^5(1 — ^) comes from the virtual 
gluon loop diagrams and the term which is suppressed by l/N^, is the quark loop contribution. Now we are ready to 
remove the collinear singularities by redefining the gluon distribution and the gluon fragmentation function as follows 



gix,^)=g^">ix)~^ 



1 asifJ.) [^ di 



I/h^Mf) 



e 2tt 



(70) 
(71) 



which is in agreement with the DGLAP equation for the gluon channel. 

Now we are ready to assemble all the rest of the finite terms into the hard factors. Let us take the finite terms 
left in the virtual contribution as an example. Using Eq. pip to perform the ^j^ integration, the finite terms are 
proportional to 



/ (fqi±(fq2±K.{qi±,q2±,q2± - k±) 



In :^ + In Ik_^^^ + "^2^ " "^i^)' 



The evaluation of the first term is trivial since it is independent of ^, qi^. Using Eqs. (|M1[55|) . the second term yieldsLl 



(27r) 



xin 



S^^\b^~x^) d\ 



p~i^'kj_-{hj_-x^) 



^ik± -r I 



{bi_ - x^f 



(73) 



Summarizing the above calculations, for the gluon channel contribution: gA — > h/ g + X, we find that the factor- 
ization formula can be explicitly written as 



^Z^p+A^h/g+X 

dyd'^p± 



/dz dx 
-^—£.xg{x,ti)Dh/g{z,n) 

d'^x±d'^y± (2) 



(2) 



(27r) 
d^x±d^y±d%± ^(2) 



Sy {x±,yi_)SY iy±,x^) 



299 ^ 27r 299 



i2.r 



S'y'ix^,b^)s'^\b^,y^)^Hg\ 



d^X^d^y^d%x_ (2) Nr.(2)., xc.(2)^ N"^a7(l) I 

—J S'y' {x±,bA_)S\^' {b±,yi)S\r' {yi^^x^)—Hl^g '> . 

The leading order results have been calculated as shown in Eq. ([S]) , from which we have 

where k± = p±/ z and r±=x±—y±. It is straightforward to show that n2qq and Hg„„ read as follows 



(74) 



(75) 



<l=^c^99(01n 






gg 

11 mfTR\ 
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NJ{l~^)e-'''^-''^ln 
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Jo 
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(76) 



(77) 



"^ The expression in Eq. I I73I I looks slightly different from the final results as shown in Eq. jTSj. Since the S- matrices are symmetrical 
among all the transverse coordinates which are all integrated over in the end, one can exchange the definition of variables xj^ <-> J/± and 
reverse the orientation of all the coordinates in Eq. I |73ll . This allows us to show that these two expressions are equivalent. 
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and 
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(78) 



Again, by choosing /i = co/r+ for the factorization scale, we can further simphfy 'H2„„ and obtain H 



2gg 



C. The quark to gluon channel q ^^ g 

This channel is relatively simpler than the q ^ q channel for two reasons: first, there is no virtual graphs; second, 
there is no rapidity divergence in the real contributions since the lower limit of the gluon longitudinal momentum is 
bounded by the hadron longitudinal momentum. It is quite straightforward to write down the cross section for this 
process by integrating out the phase space of the final state quark (/c2 , k2±) in Eg. dlip . Then, we can transform the 
cross section into momentum space and take the large iVj, limit. In the end, we obtain 



, pA^h/g+X 
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' d^ 



dz 
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k± - S,q2J 



{ki_ - qii_ - q2A 



{ki_ - ^92J 



(79) 



where we have defined ^ to be the longitudinal momentum fraction of the gluon with respect to the initial quark. 
The production of small-x gluon in pA collisions has been studied quite some time ago in Ref. [4l|. We find com- 
plete agreement between our calculation and the partonic results in Eqs (56-58) of Ref. [41| if we remove the gluon 
fragmentation function and the quark distribution, and take the limit ^ — > with dy = -^. 



Following the same procedure, we remove the collinear singularities as follows 

5(x,,).,(0)(.)-i^^('^) f''^' 
e 



j'^CFV,,{Oq 



(80) 
(81) 



where we renormalize the gluon distribution and quark fragmentation function in this off-diagonal channel. Here we 
have defined Vgq{S,) = Ml + (1 — C)^] and substituted -^ by Cp since they are equal in the large Nc Hmit. 
Therefore, we find the factorized cross section in this channel can be written as 
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(82) 



By defining W{ki^,k2±) = e-*'=i^-(^^-!'^)-''=2i-(2'^-^^), we find 
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(83) 



We can also choose fj, = co/r± for the factorization scale which yields H 



(1,1) 

2S9 



H 



(1,2) 

2sg 
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D. The gluon channel g 



To complete the calculation for all the channels, we should compute the g — > q q, although it is suppressed by a 
factor of -j^. For the gluon channel g — ?► qq, we can start from Eq. (88) of Ref. [S^l, which allows us to obtain 
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Following the above procedure, we remove the collinear singularities as follows 
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(84) 

(85) 
(86) 



where we renormalize the quark distribution and gluon fragmentation function in this off-diagonal channel. Here 



n.(0= (i-o' + e' 



In the end, the factorization formula for the cross section is 
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IV. CONCLUSION 



In summary, we have calculated the NLO correction to inclusive hadron production in pA collisions in the small- a:; 
saturation formalism. The collinear divergences are shown to be factorized into the splittings of the parton distribution 
from the incoming nucleon and the fragmentation function for the final state hadron. As we have shown above, the 
renormalization of the parton distributions of the proton and fragmentation functions follow the well-known DGLAP 
equation 
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(89) 



(90) 



In the dimensional regularization, the most common convention for the gluon spin average is to use „,._ , = -^{1 + e + ■ ■ ■). The term 
which is proportional to e can combine with the - pole terms and give a finite contribution as seen in the second term in the square 



brackets in T-i 



(1,1) 

299 



and H 



(1,2) 
299 ■ 



19 



respectively. The rapidity divergence at one-loop order is factorized into the BK evolution in either fundamental 
representation or adjoint representation for the dipole gluon distribution of the nucleus. The hard coefficients are 
calculated up to one-loop order without taking the large Nc limit for the quark q —> q channel. For some technical 
reasons, especially avoiding the sextupoles, as we have explained during the derivation, we take the large Nc limit 
for other channels. In principle, using these hard coefficients together with the NLO parton distributions and frag- 
mentation functions as well as the NLO small- a: evolution equation[42, |43| for dipole amplitudes, one can obtain the 
complete NLO cross section of the inclusive hadron production in pA collisions in the large Nc limit. The corrections 
to this NLO order cross section are either of order a^ or suppressed by j^. As to the running coupling effects 44] in 
our hybrid factorization formalism, we have no as dependence at the leading order {ag has been absorbed into the def- 
inition of the saturation momentum), and one power of a^ at the NLO, thus we find that the one-loop approximation 
for the running coupling should be sufficient. 

We have shown that the differential cross section for inclusive hadron productions in pA collisions can be written 
in a factorization form in the coordinate space. The factorization scale dependence in the hard coefficients reflects 
the DGLAP evolutions for the quark distributions and fragmentation functions. It is interesting to note that similar 
coordinate dependence (associated with r±) has also been found in the transverse momentum resummation formalism 
derived for the Drell-Yan lepton pair production in Ref. [4^. On the other hand, the hard coefficients in our case do 
not contain double logarithms, therefore there is no need for the Sudakov resummation for forward inclusive hadron 
production in pA collisions. 

Adding all the channels together in the large Nc limit gives 
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with factorization scale chosen as /^ = co/ri_ and 
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(95) 



where all the hard factors are defined in previous section. Since now the factorization scale /i depends on r±_ , the parton 
distributions and fragmentations function should change accordingly when we integrate over all the coordinates. In 
other words, the above expression should be understood as if the parton distributions and fragmentation functions 
are written inside those coordinate integrals. 

In addition, we have also demonstrated that all the hard factors can be calculated easily in the well-known MV and 
GBW model and shown that our results agree with the coUinear factorization results in the dilute limit. 

In the above calculations, we focus on the hadron production in the forward p^ collisions, where we can safely neglect 
the transverse momentum effects from the incoming parton distributions of the nucleon. The explicit calculations at 
one- loop order in the above also support this factorization, i.e., the colhnear divergence associated with the incoming 
parton distribution from the nucleon does not contain the transverse momentum dependence. The situation may 
change if we have both small-x effects from nucleon and nucleus, such as in the mid-rapidity in pA collisions at the 
LHC, when the transverse momentum effects from the gluon distribution of nucleon become important. It is in this 
region that a naive fc_L-factorization has been derived (2|, [j] and has been widely used in the literature. It will be 
interesting to extend our calculations to this kinematics too. We leave this for a future publication. 
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